/************************************************************************
*************************************************************************
*************************************************************************

simulate_example_datasets.do

Very simple script to generate a number of simulated datasets that
can subsequently be used for estimation or for calculating summary stats.

20200123: created

		
*************************************************************************
*************************************************************************
************************************************************************/


clear
clear matrix
set more off


if "`logon'" == "1" {
	cap log close

	cap mkdir "logs"
	local tmp = string(date("`=c(current_date)'", "D M Y"),"%tdCCYYNNDD")
	cap mkdir "logs/`tmp'"
	local tmp2 = clock("`=c(current_time)'","hms")
	log using "logs/`tmp'/`tmp2'simulate_example_datasets"
}


//Load globals
do "constants/simglobals.do"


//Targets for fraction of firms constrained
global consttargets 80 70 60 50 40 30 20 10 1 0.01

local refreshcon 1

//Name of folder where we'll save the results
global savename exdata

set seed 123456789

/************************************************************************
*************************************************************************
Part 1: Helper Programs
************************************************************************
************************************************************************/

//Load data generator
do "do/programs/gendatsyn.do"


/************************************************************************
Program: product
Purpose: Computes binary interaction terms
************************************************************************/
cap program drop product
program product
	syntax varlist
	
	gettoken firstvar resvars : varlist	
		
	//Take the product
	foreach var of varlist `varlist' {
		cap gen `firstvar'X`var' = `firstvar'*`var'
	}
	
	//If we're not on the last word, keep going
	if !mi(word("`resvars'",1)) {
		//Continue recursion
		product `resvars'
	}
end




/************************************************************************
Program: truemoments
Purpose: Computes true moments (must tell it the true production function)

Syntax:
name	-	"ces" "quad" (anything else = cobb douglas)
pm		-	Parameter on M (theta)
pk		-	parameter on K
pl		-	parameter on L
mu		-	elasticity of substitution
************************************************************************/
cap program drop truemoments
program truemoments
	syntax name(name=name), pm(real) pk(real) pl(real) [  mu(real 0) sc(real 0) ]
	
	/************************
	Compute marginal products
	***********************/
	
	/****** CES ********/
	if "`name'" == "ces" {
		//For CES they need to pass mu
		
		if mi("`mu'") {
			error "Must pass mu (elasticity of substitution) for a CES production function"
		}
		
		if mi("`sc'") {
			local sc = 1-`pm'
		}
		
		/* There seems to be an error here
		//Marginal product of K and L
		foreach X of varlist K L {
			local x = lower("`X'")
			gen MP_`X' = `p`x''*(1-`pm') * Y^(  (1-`pm'-`mu')/(1-`pm')  ) * exp(omega + eps)^( `mu'/(1-`pm') ) / `X'^( 1 - `mu' )
		}
		*/
		
		tempvar I
		gen `I' = `pk'* K^`mu' +   `pl' * L^`mu'
		
		foreach X of varlist K L {
			local x = lower("`X'")
			
			
			
			//This is actually the elasticity times Y/X
			gen MP_`X' =  `sc' *`p`x'' * `X'^`mu' / `I'  * Y/`X'
		}
		
		//Marginal product of M is just the cobb-douglas expression
		gen MP_M = `pm' * Y/M
		
		di "This is CES..."
		
		
	}
	
	/****** Quadratic ********/
	else if "`name'" == "quad" {
		//Marginal product for K and L is just param * Y/X
		foreach X of varlist K L {
			local x = lower("`X'")
			gen MP_`X' = `p`x'' * Y/`X'
		}
		
		//Marginal product of M is slightly different
		gen MP_M = exp(omega + eps) * K^`pk' * L^`pl' * (`pm' - 2*M)
	}
	
	/****** Cobb-Douglas ********/
	
	else {
		//Marginal product is just param * Y/X
		foreach X of varlist K L M {
			local x = lower("`X'")
			gen MP_`X' = `p`x'' * Y/`X'
		}	
	}
	
	
	
	
	/************************
	Moments: Elasticities and Sdev of marginal products
	***********************/
	
	//Elasticity = MP_X * X/Y
	foreach X of varlist K L M {
		local x = lower("`X'")
		
		//This may vary by firm
		gen elas`x'_tru 	= MP_`X' * `X'/Y 
		
		//Exclude the top and bottom 1 percent
		bysort t: cumul MP_`X', gen(MPpct)
		
		bysort t: egen sdevMP`x'_tru 	= sd(MP_`X') if MPpct > .01 & MPpct < .99
		
		//Clean up
		drop MP_`X' MPpct
	}

	
end



/************************************************************************
Program: getconstraints
Purpose: Sets a global equal to the proper values of the constraint set

Arguments:
spec	-	The specification to use

************************************************************************/
cap program drop getconstraints
program getconstraints
	args spec
	
	quietly {
	
		
		preserve
		
		local count = 0
		foreach const of numlist -1(.2)25 {
			spec`spec' 2^`const'

			collapse constrained
			
			
			
			local count = `count' + 1
			local frac`count' 	= constrained in 1
			local const`count' 	= `const'
		}

		clear
		set obs `count'
		gen frac 	= .
		gen const 	= .
		forvalues i=1/`count' {
			foreach var of varlist _all {
				replace `var' = ``var'`i'' in `i'
			}
		}
		
		

		global const`spec' ""
		local conlist "${consttargets}"
		foreach target of numlist `conlist' {
			gen tmp = abs(frac-`target'/100)
			sort tmp
			local tmp1 = const in 1
			
			global const`spec' "${const`spec'} `tmp1'"
			drop tmp
		}
		
		restore
		
	}
end






/************************************************************************
*************************************************************************
Part 2: Scenarios
************************************************************************
************************************************************************/


/************************************************************************
Fixed parameters
************************************************************************/

//List of constraints (in log base 2)
//local conlist -3(2)7

//Parameters of the production function
global pk 		= .11
global pl 		= .28
global pm 		= .67


//global quad 	= 10


//CES, almost cobb-douglas
global mu 		= .2   //Implies an elasticity of substitution equal to 1.25
global cespk 	= .064
global cespl 	= .33
global sc 		= .39

//CES, not very cobb-douglas
global mu2		= .8   //Implies an elasticity of substitution equal to 5
global cespk2	= .015
global cespl2	= .95
global sc2		= .39

////Standard deviations

//Anticipated shock
global sigeta1 	= .0925324 
global sigeta2 	= ${sigeta1} * 2

//Unanticipated shock
//global sigeps 	= 1.07
global sigeps 	= .2542193



//Number of firms
global nfirmnorm 	2613
global nfirmsmall 	1000


//Number of years
global nyearnorm 	7
global nyearshort 	4

//Persistance of productivity
global rhocal .7709392



/*********************
Production Functions
**********************/

/*
The DGP program can handle production functions of the form

X(K,L)*M^pm

OR

X(K,L)*(pm1*M - M^2)

where X(K,L) is arbitrary. This works because in both cases optimal M depends on
only X(K,L), but not K and L separately.
*/

global XKL1 "gen X_kl = K^${pk} * L^${pl}"
global XKL2 "gen X_kl = (  ${cespk}*K^${mu} + ${cespl}*L^${mu} ) ^ ( ${sc}/${mu} )"
global XKL3 "gen X_kl = (  ${cespk2}*K^${mu2} + ${cespl2}*L^${mu2} ) ^ ( ${sc2}/${mu2} )"


/*********************
Productivity processes
**********************/
global nlomproc "86.30933 - 7.081727 - 33.35803*(l.omega + 7.081727) + 4.568181*(l.omega + 7.081727)^2 - .2030311*(l.omega + 7.081727)^3 "



/************************************************************************
Specifications
************************************************************************/
local nspecs = 1

cap program drop spec1

//Baseline
program spec1
	args const fracobs
	

	
	gendatsyn  =`const' , klfunc("${XKL3}") rho(${rhocal}) n(${nfirmnorm}) t(${nyearnorm}) pm(${pm}) sigeps(${sigeps}) sigeta(${sigeta1}) 
	
	truemoments ces, pm(${pm}) pk(${cespk2}) pl(${cespl2}) mu(${mu2}) sc(${sc2})
end




/************************************************************************
*************************************************************************
Part 3: Sweep
************************************************************************
************************************************************************/
cap mkdir "dta/simulations"
cap mkdir "dta/simulations/${savename}"

if ${ncores} > 1 {
	cap mkdir "dta/simulations/${savename}/1_1"
	
	local savefile "dta/simulations/${savename}/1_1"
	
	di "*******************"
	di "I am job `jobid'"
	di "*******************" _n _n _n
}
else {
	local savefile "dta/simulations/${savename}/1_1"
	cap mkdir "`savefile'"
}



//Loop through specifications
foreach spec of numlist 1 {
	di _n "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"
	
	
	//Get the constraints for this specification
	getconstraints `spec'
	local conlist "${const`spec'}"
	local ind 1
	foreach const of numlist `conlist' {
		local fc : word `ind' of $consttargets
		
		cap mkdir "`savefile'/const`fc'"
		
		di "-----------------------------------------------------------------------------------------"
		forvalues b=1/${Bser} {
			
			di "On Specification `spec'" _col(25) "Constraint = `fc'" _col(70) "Rep = `b'"
			quietly {
				spec`spec' 2^`const' 
				order elas* sdevMP*		
				
				save "`savefile'/const`fc'/`b'.dta", replace
				
			}
			
		}
		
		local ++ind
	}
	
}
